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Abstract Sixty days of Doppler images from the Solar and Heliospheric Obsen’atory 
(SOHO) / Michelson Doppler Imager (MDI) investigation during the 1996 and 2008 so- 
lar minima have been analyzed to show that certain supergranule characteristics (size, size 
range, and horizontal velocity) exhibit fluctuations of three to five days. Cross-correlating 
parameters showed a good, positive correlation between supergranulation size and size 
range, and a moderate, negative correlation between size range and velocity. The size and 
velocity do exhibit a moderate, negative correlation, but with a small time lag (less than 12 
hours). Supergranule sizes during five days of co-temporal data from MDI and the Solar 
Dynamics Observatory (SDO) / Helioseismic Magnetic Imager (HMI) exhibit similar fluc- 
tuations with a high level of correlation between them. This verifies the solar origin of the 
fluctuations, which cannot be caused by instrumental artifacts according to these observa- 
tions. Similar fluctuations are also observed in data simulations that model the evolution of 
the MDI Doppler pattern over a 60-day period. Correlations between the supergranule size 
and size range time-series derived from the simulated data are similar to those seen in MDI 
data. A simple toy-model using cumulative, uncorrelated exponential growth and decay pat- 
terns at random emergence times produces a time-series similar to the data simulations. The 
qualitative similarities between the simulated and the observed time-series suggest that the 
fluctuations arise from stochastic processes occurring within the solar convection zone. This 
behavior, propagating to surface manifestations of supergranulation, may assist our under- 
standing of magnetic-field-line advection, evolution, and interaction. 
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1. Introduction 

Irregularities in the measured rotation rate of the Sun were noticed in the early 1900s, 
but were attributed to observer bias (see, inter alia, Plaskett and Delury, 1913; Plaskett, 
1916). Hart (1956) showed that the fluctuations were solar in origin and had scale lengths 
of ~ 26000 km. It was the use of spectroheliograms by Leighton, Noyes, and Simon (1962) 
that allowed the Doppler velocity to be measured across the disk of the Sun and showed 
that these velocities formed the cellular patterns we call supergranules (see the review by 
Rieutord and Rincon, 2010). 

Supergranulation cells have lifetimes of one to two days and are typically ~ 30 Mm 
across, but the range of sizes can be quite large (Meunier et al., 2007; Rieutord et al., 2010). 
Table 1 summarizes the findings of Williams and Pesnell (2011a). It was found that dur- 
ing both the most recent solar minima (that concluded Solar Cycles 22 and 23, respectively) 
supergranulation cells were on the order of 35 Mm across. Using the upper and lower bound- 
aries of the FWHM of the supergranule feature in the cell-size distribution spectrum, a size 
range between 19-65 Mm was estimated. Previously, Hathaway et al (2000) produced size 
distributions of Doppler patterns within high-resolution Michelson Doppler Imager (MDI) 
images that exhibit a crossover between the supergranule and granule components. This sug- 
gests that the lower limit of supergranulation cell sizes may extend to below 10 Mm, similar 
in size to the largest granules. 

Supergranules are well-known for their low radial velocities and strong divergent hor- 
izontal surface-flows. By isolating the spectral power at supergranule wavenumbers, thus 
removing the velocity contributions from granulation, Williams and Pesnell (2011a) found 
the radial and horizontal velocities to be around 8 ms -1 and 140 ms -1 , respectively (Ta- 
ble 1). Although it seems that little heat is convected by these cells, reflected in the observed 
weak center-to-edge temperature gradient of photospheric supergranulation cells (Meunier, 
Roudier, and Rieutord, 2008; Goldbaum et al., 2009), they certainly play a role in dragging 
magnetic-field lines to the surface and advecting them to cell boundaries. 

While the aforementioned studies of the long-term, inter- year variation of supergranula- 
tion characteristics exist, studies at much shorter timescales are rarer. We have extended our 
studies (Williams and Pesnell, 2011a) of supergranulation characteristics (such as sizes, size 
ranges, and velocities) at solar minimum to provide a better error- and statistical analysis of 
our data, and in doing so have observed fluctuations in these characteristics with periods of 
three to five days (Williams and Pesnell, 2011b). This article describes the observation and 
subsequent analysis of the time-series of these properties and any correlations among them, 
and also determines whether the sources of the fluctuations are instrumental or solar. 


2. Data Preparation 

The Michelson Doppler Imager (MDI: Scherrer et al., 1995), launched onboard the Solar 
and Heliospheric Observatory (SOHO) in 1996, has performed a yearly dynamics mn of 
approximately 60 days during which full-disk 1 024 x 1024 Dopplergrams were transmitted 
to the ground. These 2" pixel -1 resolution image data contain line-of-sight photospheric 
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Table 1 Comparison of 
time-averaged supergranule 

Parameter 

Symbol 

Units 

1996 Mean 

2008 Mean 

parameters between the minima 
of Solar Cycles 22/23 and 23/24. 

Spectral peak wavenumber 

^peak 


122 ± 1 

125 ± 1 


Spectral peak FWHM 

Al 


155 ±3 

164 ±4 


Supergranule diameter 

X 

Mm 

35.9 ±0.3 

35.0 ±0.3 


Radial velocity 

V T 

ms -1 

7.6 ± 1.4 

7.5 ±1.6 


Horizontal velocity 

v h 

ms -1 

139 ±1 

141 ± 1 


1/e Lifetime 

T l/e 

h 

17.8 ±0.5 

17.7 ±0.4 


Number of supergranules 

N SG 


6005 

6318 


velocities obtained from the Ni I 67 67 A line at a one-minute cadence. These images have 
been analyzed to produce an array of results that include helioseismic studies of the internal 
differential rotation profile (Schou et al., 1998) and determinations of subsurface flows using 
methods such as ring-diagram and time-distance analyses (Hindman et al., 2004). 

For our studies of supergranulation, these images require processing to remove the 
Doppler signals caused by the five-minute p-mode oscillations that have spatial scales 
similar to supergranules. The removal process (Hathaway, 1988b) involves averaging 31 
Gaussian-weighted Dopplergrams that were derotated to the rotational phase of the central 
image. The resulting averages are then sampled every 15 minutes, producing 96 Doppler 
images per day. 

The SOHO spacecraft orbits around the Li Lagrangian point situated at the point of 
gravitational balance between the Sun and the Earth. Along with the yearly orbital motion 
around the Sun, the velocity signal due to the spacecraft’s motion around Li needs to be 
removed from the Dopplergrams. Furthermore, the gravitational redshift velocity component 
is also subtracted. 

Finally, the axisymmetric flows, such as differential rotation, meridional flow and convec- 
tive blueshift, are removed. The method, described by Hathaway (1992), involves fitting the 
flow profiles by calculating the coefficients of a spherical-harmonic series that are summed 
to give the velocity profile of each flow and then subtracting that profile from the image. 
After also removing any instrumental artifacts, the resulting image contains velocity signals 
due only to photospheric manifestations of solar convection. An additional step during this 
process involves remapping the image to heliographic coordinates (Thompson, 2002), linear 
in both latitude and longitude. 


3. Data Analyses 

We applied analysis methods similar to those used by Williams and Pesnell (201 la) to cal- 
culate from each processed Doppler image the globally averaged supergranule spectral peak 
wavenumber [£ peak ] and the FWHM | Af ], from which the supergranule size and size range 
can be derived, as well as the horizontal- velocity component [ F h ] . Data between 25 May and 
23 July 1996 (MDI day 1240- 1299) and 3 March and 1 May 2008 (MDI day 5540- 5599) 
were analyzed. 

To determine £ pea k and Al, a modified Forentzian was used to fit the supergranule feature 
of the power spectrum derived from each heliographic Doppler map. Spectral amplitudes 
[A" ! ] in spherical harmonic degree \ l \ and order [m] were derived by projecting the maps 
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Spherical Harmonic Wavenumber 



sin p 


Figure 1 Analysis plots of a single MDI Dopplergram from 27 April 2008, 00:00 UT. (Left) A power 
spectrum derived from a Doppler image (black dots) is fitted with a modified Lorentzian (solid red line). 
From this, f pea k and At (signified by the FWHM, dashed red line) are determined. The full function is 
plotted for context (dotted blue line). (Right) Average values of the squared velocities of pixels within each 
annulus at an angular distance [p] from disk center are plotted against sin 2 p. The horizontal velocity [ Vj-,] is 
determined from the fit. 


onto the spherical harmonics (Hathaway et al., 2000). The power spectrum in i was then 
produced by summing over m : 


i 

no= EK 

m=—i 


(i) 


The supergranule feature, as seen in Figure 1 (left) peaking at around t = 120, is fitted 
between 50 < l < 250 with a modified Lorentzian, 


m = 


IA 

[(f-f 0 ) 2 + T 2 J 


(2) 


from which l pea k is derived from the peak of the fitting function and At from its FWHM. 

The horizontal velocities [ Vfi] were determined from each heliocentric Dopplergram by 
dividing the image into a set of concentric annuli of an angular distance [p] from disk 
center and then determining the average of the squared velocities at each pixel contained 
within (Hathaway et al, 2002). Plotting the average value within each annulus against sin 2 p 
produces a straight line (Figure 1 , right) from which Vh can be calculated from the slope and 
intercept values using 


V 2 (p) = V} + ( V h 2 - V}) sin 2 p. 


( 3 ) 


However, to properly determine Vh throughout the time-series, Doppler image pixels 
within regions of high magnetic activity must be removed because spurious velocity values 
are measured in these regions (Liu and Norton, 2001). We used contemporaneous 96-minute 
cadence MDI magnetograms to mask out Dopplergram pixels, correlating to magnetogram 
pixels with line-of-sight magnetic fields above 25 Gauss. We then used morphological open 
and close algorithms to mask out pixels over a large area to make sure that any spurious 
velocity values not matching the original mask were also masked out. We then applied the 
horizontal- velocity calculations to each individual Dopplergram. 

The resulting data consist of three 15-minute-cadence time-series for each of the quanti- 
ties analyzed. The mean and standard deviation for each time-series is calculated and shown 
in Table 2. 
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Table 2 Comparison of supergranule parameters between the solar minima of Cycles 22/23 and 23/24 de- 
termined by averaging over the 15-minute-cadence time-series. While the mean values are similar to those in 
Table 1, the error estimates are now represented by the standard deviations of the datasets. 


Parameter 

Symbol 

Units 

1996 Mean 

2008 Mean 

Spectral peak wavenumber 

^peak 


122 ± 1 

125 ±2 

Spectral peak FWHM 

M 


156 ±7 

164 ±9 

Horizontal velocity 

v h 

ms -1 

139 ±3 

141 ±3 



Days 

Figure 2 A time-series of fp ea k calculated from the spectra derived from the set of 1996 Dopplergrams. 
Each 15-minute £p ea k value is represented by a black dot. The dashed line corresponds to a mean value 
Zpeak = 121.7, the dot-dashed and dotted lines are the 1 a and 2 a values, respectively, where a = 1.5. The 
blue line is due to smoothing the time-series with a 33-point-wide Savitzky-Golay filter that retains much of 
the variation of the original time-series. 


An example time-series for l pea k is shown in Figure 2 along with the values for the mean 
[fpeak = 121.7] and standard deviation [cr = 1 .5]. The time-series exhibits fluctuations with a 
period of around four days, which were investigated by means of frequency analysis. Firstly, 
the data were smoothed using a fourth-degree zeroth-order 33-point-wide Savitzky-Golay 
filter (overplotted in Figure 2). The mean was subtracted and low-frequency trends were re- 
moved by subtracting a version of the original time-series that was smoothed with a five-day 
boxcar function. These trends were not investigated because the current work focuses on the 
four-day periodicities. After padding the ~ 6000 point data-series to 8192 points to reduce 
aliasing and normalizing to a standard deviation of one, a fast Fourier transform (FFT) was 
applied, resulting in the frequency spectra shown in Figure 3. The window function illustrat- 
ing the resolution of the frequency analysis is shown as a dashed black line. The frequency 
distribution in Figure 3a relates to the time-series of Figure 2. Dominant peaks are seen at 
around 0.25 and 0.35 cycles per day. Peaks at similar frequencies (and thus similar periods) 
are also seen in the 2008 data (Figure 3b). Fourier analyses were also applied to the 1996 
and 2008 time-series for At and F h with si mi lar results (Figure 3c -f). These show that 
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Figure 3 All of the filtered time-series are de- trended, have their means removed, are padded to 8192 points, 
and normalized to a standard deviation equal to unity before applying an FFT. The resulting 1996 (a, c, e) 
and 2008 (b, d, f) power spectra (solid blue line) are plotted for the supergranule characteristics f pea k, A£, 
and Vh, respectively. A plot of the window function (dashed black line) is added, representing the resolution 
of the frequency analysis and peaking at an arbitrary frequency for clarity. 


there seems to be an underlying time-dependent fluctuation for each of these characteristics 
of around 0.2 - 0.4 cycles per day, relating to periods of around three to five days. 


4. Inter-Characteristic Comparisons 

With the observed fluctuations in the time-variation of each of the supergranule character- 
istics, it is worth investigating any correlation between them. The Pearson correlation co- 
efficient, which parameterizes their mutual correlation, was calculated for each time-series 
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Table 3 Comparison of 
time-averaged supergranule 

Quantity 

Characteristics 

1996 

2008 

parameters between the solar 





minima of Cycles 22/23 and 

Pearson correlation coefficient 


0.41 

0.48 

23/24. 


peak ’ Hi 

-0.03 

-0.10 



R At,V h 

-0.21 

-0.30 


Hurst exponent 

^peak 

0.76 

0.80 



A l 

0.82 

0.82 



Vh 

0.76 

0.80 


pair. For two arbitrary time-series [x] and [y] consisting of N x and N y values [x,] and | y, ] , 
the Pearson coefficient for the correlated pair [/?* y ] is calculated from 



where x, y, a x and a y are the means and standard deviations of the x and y time-series, 
respectively. Note that the time-series components for this calculation include points for 
which values exist at the same time in both series. If no point exists at a certain time in 
either series, the points were removed from the calculation. 

The Pearson coefficients for each parameter pair for both datasets are listed in Table 3. 
They show that the spectral peak wavenumber (inversely proportional to the globally aver- 
aged supergranule size) and the spectral FWHM, the supergranule size-range, are moder- 
ately correlated at /?£ pcak , Ae = 0.41 for 1996 and Ri peak , m = 0.48 for 2008, while the peak 
wavenumber and horizontal velocity are poorly correlated, suggesting no relation between 
these two quantities with Rt peak ,v h = —0.03 for 1996 and ^ peak ,v h = —0.10 for 2008. How- 
ever, the supergranule size range and horizontal velocity show a weak correlation, with 
R a i v h = —0.21 for 1996 and R Af Vh = —0.30 for 2008. Figure 4 shows a set of scatter plots 
illustrating the time-series pair correlations. The mean value and lcr and 2cr values for each 
time-series are also shown. 

To consider the possibility of phase shifts between the correlated pairs, the central 30 
days of the first time-series were selected, then the Pearson coefficient calculated over a 
series of lag-times (between ±15 days) with the second time-series. Figure 5 shows the 
change in the Pearson coefficient over the range of lag-times for each parameter pair. The 
plots show that the correlation between f peak and At is strongest at zero time-lag, showing 
that there is no lag between the quantities and suggesting one quantity does not influence 
another. Both these quantities seem to be roughly anti-correlated with V h and exhibiting 
an element of lag. With t pea k being anti-correlated to V h and inversely proportional to the 
supergranule size, the results suggest a direct correlation between supergranule size and 
horizontal velocity. This matches previous findings (Meunier et al., 2007). 

These resulting inter-correlations could be caused by the size parameters being derived 
from the same analysis process while a different process was used to derive the velocities. 
They also show fluctuations in the Pearson coefficients with periods similar to the fluctua- 
tions of the time-series themselves. This is very noticeable in the 2008 data, where fluctua- 
tions in the correlation coefficient of around 3.5 days are apparent. 
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Figure 4 Scatter plots illustrating the correlations between the 60-day time-series of supergranule charac- 
teristics for 1996 (left column) and 2008 (right column). The plots correspond to the Pearson coefficients 
Rg peak, At! (top), ^£ peak ,V h (middle) and RAl,V h (bottom), respectively. The dashed, dot-dashed, and dotted 
lines relate to respective the mean, lc, and 2 a values for each time-series. 


5. Additional Analyses 

To quantify any stochastic element of the fluctuations, and thus a level of predictability, the 
Hurst exponent (Mandelbrot and Wallis, 1969; Pesnell, 2012) was determined for each of 
the six time-series. A value of around 0.5 relates to a pure random walk, while a value to- 
ward unity defines fully predictable time-series. Table 3 lists the calculated Hurst exponents, 
indicating some predictability within the time-series. 

To investigate whether the observed fluctuations are solar in origin or are due to instru- 
mental artifacts of MDI, such as defocusing (Korzennik, Rabello-Soares, and Schou, 2004; 
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Figure 5 Dependence of correlation values on lag-times between two 30-day selections from time-series 
pairs from 1996 and 2008. Note that the zero time-lag correlation values are different from those of Figure 4, 
which were calculated with the full 60-day time-series. 


Williams and Pesnell, 2011a), contemporaneous Dopplergrams taken between 6 and 10 July 
2010 by both MDI and HMI were analyzed and compared. Plotting the time-series of f pea k 
obtained from both datasets (Figure 6, left) illustrates very similar trending. Note that the 
discrepancy in the mean values is caused by the different resolution of the instruments 
(Williams et al., 2013). Correlating the data sets (Figure 6, right) results in a Pearson coef- 
ficient of 0.78, indicating a very good correlation. This highlights that the same fluctuating 
phenomenon was observed in both instruments and allows us to conclude that the fluctua- 
tions have a solar origin. 
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Supergranule Peak Values 



Figure 6 (Left) Time-series of £p ea k derived from contemporaneous MDI (red) and HMI (blue) Doppler- 
grams covering the time range 6-10 July 2010. (Right) Scatter plot representing the correlation between the 
MDI and HMI time-series. The Pearson correlation coefficient is 0.78. 


6. Computer Simulations 

The fluctuations were reproduced in simulated data with the aim of shedding more light on 
their origins. Simulations that produce realistic Doppler images (Hathaway, 1988a; Hath- 
away et al., 2010) were used to produce 60 days of images at a 15-minute cadence. The 
images were produced from a synthetic input spectrum containing both supergranule and 
granule components. These provide the spectral amplitudes for the orthogonal velocity ele- 
ments that are constructed with the use of the spherical harmonics, which in turn were com- 
bined and projected onto the line of sight. The simulation algorithm advects the Doppler 
pattern according to differential rotation and meridional circulation profiles, the former ad- 
verting the cells at rates relative to their size, i.e. dependent on the depths at which they are 
anchored. Their physical size also determines the rate at which the cells evolve. 

The simulated dataset was then processed to remove the axisymmetric-flow signals, 
mapped to heliographic coordinates and projected onto the spherical harmonics from which 
the spectral amplitudes were extracted. From these, we produced power spectra and calcu- 
lated the supergranule-feature peak wavenumber and widths. Figure 7 (top-left) shows the 
temporal distribution of the peak wavenumbers, showing fluctuations that seem to have a 
higher frequency than those seen in the real data. This is highlighted in the spectrum (Figure 

7, top-right) of the time-series, where a prominent peak is observed at 0.5 cycles per day, 
i.e. a periodicity of two days. This is shorter than that observed in the MDI results, possi- 
bly due to the input function that determines the evolution of the convective Doppler pattern 
that overcompensates for smaller, shorter- lived features. These features tend to dominate the 
time-series, which is reflected in the power spectrum. This overcompensation is seen when 
comparing the decay of the synthetic Doppler signal with that from MDI, in particular at 
low time-lags. This is also seen at high time-lags, but by this time the longer-lived pattern 
is a much weaker component of the pattern as a whole, and so the short-lived signal is more 
dominant within the time-series and the related spectra. While this difference in periodicity 
between the observed and simulated time-series could be explained by the decay rate of the 
Doppler signal (Figure 7, bottom- left), which is an element of the simulation to be analyzed 
in future work, that the fluctuations are observed in simulated data with no hydrodynamical 
influence whatsoever suggests that they arise from the evolution of the Doppler signal itself. 
A similar trend is seen within a time-series of the FWHM values. The subsequently calcu- 
lated value of the correlation coefficient relating the two time-series (Figure 7, bottom-right) 
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Pearsons product moment coefficient = 0.43 



Figure 7 (Top-left) The peak wavenumbers for the supergranule feature in the spatial spectra derived from 
synthetic Dopplergrams (black dots) are smoothed with the same Savitzky-Golay filter as used for the MDI 
data. (Top-right) The frequency spectrum, similar to those in Figure 3, shows that, while a high back- 
ground exists, there is a dominant peak of power at frequencies relating to a fluctuation period of around 
two days. (Bottom-left) Comparison of the time-evolution of correlation coefficients tracking supergranule 
patterns within synthetic (blue) and real (red) Doppler data. (Bottom-right) The cross-correlation of the peak 
wavenumber and the FWHM of the supergranule feature time-series give similar results to those derived from 
the MDI data. 


is found to be R^ake a »fake = 0.43, where the respective mean and standard deviations for 

peak ' 

these two time-series are = 123.9 ±1.4 and Af FAKE = 167.8 ± 8.2. 

A toy model was constructed to simulate cumulative, uncorrelated exponential growth 
and decay patterns at random times. Performed over a number of iterations, a pattern 
emerges that is qualitatively similar to that of the MDI time-series (Figure 8). This simple 
model suggests that the fluctuations observed in the MDI data can be described by stochastic 
processes underlying the evolution of uncorrelated individual supergranules that comprise 
the global Doppler pattern. 


7. Conclusions and Future Work 

Time-series of globally averaged values for supergranule size, size-range, and horizontal 
velocity calculated from individual 15-minute-cadence MDI Doppler maps covering a time 
period of 60 days exhibit fluctuations with periods of three to five days. We found that the 
time-series for supergranule sizes are consistent between co-temporal MDI and HMI data, 
meaning that the observed fluctuations are solar in origin. While the source of these fluctu- 


<£) Springer 


1112 


P.E. Williams, W.D. Pesnell 



Time (days) 

Figure 8 A time-series produced by the toy model of stochastic growth and decay of supergranulation has a 
similar appearance to that derived from MDI data. 


ations is unclear, data simulations show that stochastic, uncorrelated variations seem to be 
qualitatively similar to the observational data. 

Future studies of the time-evolution of the supergranulation pattern, including its growth 
and decay, may shed light on any influences underlying the time-series characteristics. This 
can include the adjustment of the input function that determines the pattern evolution within 
the data simulations to produce decay rates that more closely follow those seen within the 
real data. 

With supergranulation being visible in other, non-Doppler datasets, such as Ca II K im- 
ages, it is of interest to study whether these fluctuation patterns manifest in such data. McIn- 
tosh et al. (2011) highlighted variations in supergranulation properties within Ca II K images 
over long timescales at large time-steps, so it will be of interest to analyze these properties 
over shorter intervals with a more rapid cadence. 
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